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order  Runge-Kutta  algorithm  incorporating  an  embedded  fourth-order  result.  This 
latter  alternative  includes  automatic  time-step  modification  and  guarantees  a 
prescribed  level  of  accuracy.  Several  utility  routines  are  provided  in  support 
of  the  method  of  characteristics.  There  is  a  two-dimensional  spline  calcula¬ 
tion  procedure  for  the  analytic  description  of  flow-field  parameters,  steady  po¬ 
tential,  and  potential  gradient.  Spline  interpolation  subroutines  enable  the 
user  to  incorporate  function  and  gradient  evaluations  directly  into  computer 
program  code.  There  is  also  a  routine  for  calculating  average  linear  velocity 
for  those  flow  situations  where  Darcy's  law  is  appropriate.  A  plotting  routine 
features  an  option  for  producing  flow-line  map  overlays  to  a  user-supplied 
scale. 
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A  PROCEDURE  FOR  CALCULATING  GROUNDWATER  FLOW  LINES 

by 

Charles  J.  Daly 

Any  continuous  mathematical  function  of  spatial  coordinates  x  and  time  t 
can  constitute  a  field.  Fields  may  be  scalar,  vector,  or  tensor,  e.g.  hy- 

-►  -y  -y 

draullc  head  h(x,t),  velocity  v(x,t),  and  stress  T(x,t)  respectively.  The 
scalar  field  is  the  most  elementary  type.  Components  of  vector  and  tensor 
fields  are  themselves  scalar  fields. 

The  principles  of  classical  physics  generally  suffice  to  describe  the 
motion  of  a  group  of  individual  rigid  objects.  For  example,  the  velocity  of 
the  center  of  mass  of  the  1th  object  in  a  group  can  be  specified  as  a  func¬ 
tion  of  time: 

V.  (t)  =  (v  (t),v  (t),v  (t)).  (1) 

1  Xi 

In  considering  fluid  flow,  however,  the  same  principles  of  classical 
physics  can  not  be  applied  to  the  uncountable  number  of  molecules  contained 
in  even  a  small  volume  of  fluid.  Instead,  the  field  approach  is  employed, 
and  the  velocity  of  fluid  particles  is  specified  by  the  velocity  vector 
field: 


V  (x,t)  =  (v^(x,t),v^(x,t),v^(x,t)),  (2) 

that  is,  a  fluid  particle  located  at  point  Xq  at  time  tg  has  velocity 
V  (Xo,to). 

The  significance  of  the  velocity  vector  field  can  be  appreciated  from 
two  perspectives.  The  first,  called  the  Eulerian  viewpoint,  focuses  atten- 
tlon  on  specific  spatial  locations  Xp  and  observes  the  velocity  of  parti¬ 
cles  moving  past  those  locations.  The  velocity  field  is  revealed  by  the  set 
of  observations: 

V  (Xp,t)  p  =  1 ,2,3,.. . 


(3) 


The  second  perspective,  called  the  Lagrangian  viewpoint,  focuses  on  a 
number  of  specific  particles  moving  with  the  flow.  Over  time  each  particle 
follows  a  trajectory,  or  flow  line,  defined  by  the  locus  of  points: 

Xj^(t)  =  (x^(t)  ,yj^(t)  ,Zj^(t))  k.  =  1,2,3,...  (4) 

Monitoring  the  velocity  of  each  moving  particle  would  give  the  set  of  data 
again  revealing  the  velocity  field: 

V  (Xj^(t),t)  k=  1,2,3,...  (5) 

The  intent  of  this  report  is  to  describe  a  methodology  for  determining 
the  flow  lines  of  particles  carried  along  as  part  of  a  moving  fluid.  It  is 
natural  then  to  adopt  the  Lagrangian  viewpoint  of  the  velocity  field. 

OVERVIEW  OF  MODELING 

Figure  1  is  a  schematic  of  the  flow-line  calculation  procedure  as  ap¬ 
plied  to  the  case  of  steady  ground-water  flow.  The  system  is  assumed  to  be 
governed  by  Darcy’s  law  and  the  mass  conservation  principle.  The  main  compo¬ 
nents  are; 

1)  the  spline  generator  and  function  subroutines, 

2)  the  average  linear  velocity  calculation  subroutine, 

3)  the  flow-line  determination  routines,  and 

4)  the  flow-line  plotting  routine. 

Figure  2  illustrates  the  results  of  applying  the  procedure  in  Figure  1. 
The  composite  map,  in  this  case  a  flow  net,  could  be  used  to  predict  the  tra¬ 
jectory  of  ground-water  contamination  and  the  propagation  rate  of  a  contami¬ 
nation  front.  (In  practice  one  must  also  consider  the  importance  of  sorption 
and  mechanical  dispersion  on  the  transport  (Daly  1983).) 

In  the  sections  that  follow,  each  component  in  the  complete  methodology 
of  Figure  1  is  described,  documented,  and  then  demonstrated. 

SPLINE  GENERATING  PACKAGE 
Description 

Consider  the  two-dimensional  domain,  shown  in  Figure  3,  defined  by 
X  £  xn  and  Yi  ^  y  Ym*  Let  D  denote  the  domain  and  3D  its  boundary. 

Values  of  a  function  f(x,y),  whose  analytic  form  is  unknown,  are  given  at  a 


Figure  1.  Schematic  of  modeling  procedure  for  two-dimensional  steady  ground- 
water  flow. 
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Figure  3.  Spline  function  domain. 


number  of  points  in  D  and  on  9d.  These  points  are  arranged  so  that  one  point 
is  located  at  the  intersection  of  every  vertical  (x  =  constant)  and  every 
horizontal  line  (y  =  constant).  Each  point  arranged  in  this  fashion  is 
called  a  knot.  (Regular  placement  of  the  knots,  as  in  Figure  3,  is  not 
essential  to  the  derivation  of  a  spline  interpolate.  This  assumption  does 
lead,  however,  to  the  relatively  simple  calculation  procedure  described 
la  ter. ) 

Let  the  approximating  expression  for  the  unknown  function  f(x,y)  be 
designated  by  S(x,y).  The  particular  form  of  S(x,y)  depends  on  the  set  of 
functions  from  which  it  is  to  be  constructed.  A  variety  of  so-called  basis 
functions  are  available  for  this  purpose.  In  cases  where  the  smoothness  of 
f(x,y)  can  be  estimated,  some  justification  can  be  made  for  the  selection 
of  particular  basis  functions.  If  f(x,y)  is  believed  to  vary  smoothly,  it 
would  appear  to  be  more  advantageous  to  construct  S(x,y)  from  functions 
exhibiting  this  property  (Rivlin  1969,  Schultz  1973). 

A  spline  is  a  composite  function.  For  domain  D  of  Figure  3,  a  bi¬ 
cubic  spline  Interpolate  S(x,y)  will  be  composed  of  (N-1)  x  (M-1)  bicubic 
polynomials.  (The  polynomials  are  bicubic  since  they  contain  powers  of 
both  X  and  y  up  to  order  three.)  Each  individual  polynomial  is  valid  only 
in  and  on  the  boundary  of  a  particular  subrectangle  of  D.  Thus  there  is  a 
one-to-one  correspondence  between  the  polynomials  and  the  subrectangles. 


Figure  4.  An  arbitrary  subrectangle  of  the  domain  C 


Each  of  the  constituent  bicubic  polynomials  has  the  same  s _  ral 
formula.  This  formula  can  be  Illustrated  by  writing  the  equation  the 
polynomial  associated  with  the  isolated  subrectangle  of  Figure  4.  In  total, 
(N-1)  X  (M-1)  of  these  equations  must  be  written  to  define  the  entire 
spline.  For  the  subrectangle  of  Figure  4,  the  polynomial  is: 

s(x,y)  =  f (x^,yj)C^(x)Cj(y)  +  f )C^(x)C^^j (y) 

+  f(x^^^,y.)C^^^(x)C.(y)  +  f (x)C.^^(y) 

+  11  (x^,yj)D^(x)C^(y)  +11  (x^,yj+l)D^(x)C.^^(y) 

^  (^+l.yj)Di+i(x)C.(y)  +  If  (Xi+i.yj+i)Di^l(x)C.^j(y) 

(x^.yj)C^(x)D.(y)  +|f  (^.yj+i)Ci(x)D.^^(y) 

If  (x.^^,y.)C^^j(x)D.(y)  +  ||  (x  ,y  (x)D  (y ) 

(x^,yj)D^(x)Dj(y)  +^  (x^,y.^^)D.(x)D.^j(y) 

2  2 

(x^^^,y.)D^^^(x)D.(y)  +  ^  (x^^l  .y^  +  j  )D  (x)D  (y) .  (6) 

The  interpolating  functions  Cj,  C^+j ,  Dj,  and  are  defined  by 

the  following  generalized  formulas: 
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Evaluation  of  the  interpolating  functions  and  their  derivatives  at  the 
knots  yields  the  important  results: 

3D  (0  1  3d  [0  J 
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The  last  four  equations  can  be  used  to  show  that  the  function  s(x,y) 

2 

has  the  following  properties:  the  values  s,  9s/ 9s/ 9y,  and  9  s/9x9y  are 

2 

equal  to  the  corresponding  values  f,  0f/0x,  0f/0y,  and  0  f/0x0y  at  each  of 
the  knots  (x^.yj),  (x^,yj+i),  (xi+i.yj),  and  (xi+i.yj+i). 

The  first  four  terms  in  each  bicubic  polynomial,  represented  by  eq  6, 
are  known  since  f(x,y)  is  given  at  every  knot.  The  spline  determination 
problem  is  then  to  estimate  the  remaining  12  terms  in  the  equation  of  each 
bicubic  polynomial.  These  unknowns  are  the  values  of  0f/0x,  0f/0y,  and 
02f/0x0y  at  the  four  corners  of  every  sub rec tangle.  They  are  found  by 
observing  continuity  and  differentiability  conditions  at  the  knots.  Since 
the  joint  edges  of  two  adjacent  sub rec tangles  share  common  knots,  the  poly¬ 
nomials  in  each  of  these  adjacent  elements  share  common  unknowns.  Tliis 
fact  suggests  the  solution  of  simultaneous  equations  for  the  unknowns,  and 
this  is  exactly  the  method  used. 


The  algebra  involved  in  obtaining  the  system  of  equations  is  quite 
cumbersome,  although  the  procedure  is  relatively  simple  to  explain.  The 
procedure  consists  of  four  steps. 


Step  1 

The  equation  of  the  bicubic  polynomial  in  each  subrectangle  is  differ¬ 
entiated  twice  with  respect  to  x.  This  yields  a  set  of  equations  for 


2  2 

9  s(x,y)/9x  ,  one  equation  for  each  subrectangle  of  domain  D. 


2  2 

A  continuity  requirement  is  imposed  on  the  determination  of  9  s/ 9x  at 


each  knot  in  D.  Since  a  given  knot  usually  belongs  to  more  than  one  sub¬ 
rectangle,  there  is  in  general  more  than  one  equation  for  the  determination 


2  2 

of  9  s/ 9x  at  a  knot.  Continuity  requires  that  all  of  these  determinations 


be  equal.  Observing  the  continuity  conditions  along  a  single  horizontal  y 
=  yj  yields  a  complete  set  of  equations  for  3f/9x  at  each  interior  knot 
on  that  line.  (The  values  of  9f/9x  must  be  given  at  the  end  points  of  that 
line.)  Each  horizontal  line  is  considered  in  succession  until  9f/9x  has 
been  determined  at  every  knot  in  D. 

Let  subscripts  i  and  J  be  used  to  refer  to  the  point  (xj.yj);  for 
example,  fj^+j  j  =  f(xi4.j,yj).  The  systems  of  equations  solved  for 
the  values  9f/9x  along  the  horizontals  can  then  be  written: 
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(16) 


S  tep  2 

A  procedure  similar  to  that  used  in  Step  1  is  used  to  determine  9f/3y 
at  every  knot  in  D.  First,  the  equation  of  each  bicubic  polynomial  is 
differentiated  twice  with  respect  to  y,  yielding  a  set  of  equations  for 
9  ^s(x  ,  y)  /  9y .  The  separate  determinations  of  9'^s/''y'^  at  a  knot  must  be 
equal.  Observing  this  requirement  along  a  single  vertical  line  x  =  x^ 
gives  a  complete  set  of  equations  for  M/^'y  along  tliat  lino. 

Given  the  values  of  9f/3y  at  the  end  points  of  each  vertical  line 
enables  the  solution  of  each  set  of  simultaneous  equations;  thus  the  values 
9f/9y  at  every  knot  are  obtained.  Tlie  systems  of  equations  can  he  con¬ 
structed  from: 


where  l^i^N,  l£j^  M-2 ,  and  Ayj  =  vj+j  -  y j . 

S  tep  3 

The  equations  of  the  bicubic  polynomials  that  are  valid  along  the 

vertical  boundaries  of  D  are  differentiated  to  produce  equations  for 
3  2 

3  s(x,y)/9x9y  .  Along  each  of  these  vertical  boundaries  the  value  of 
3  2 

9  s/9x9y  must  be  continuous  at  the  knots.  This  results  in  two  systems  of 

2 

equations,  one  for  each  boundary.  When  the  values  9  f/9x9y  are  given  at 
the  four  corners  of  D,  these  two  systems  can  be  solved  for  all  other  values 
of  9  f/9x9y  on  the  vertical  boundaries. 

The  equations  are  of  the  form: 

9^f  9^f  9^f 

2(Ay^_^,  +  AyJ  4-  Ay,  '  = 


'j+1  9x9y  ^j+] 
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9x9y 

+  Ay. 
J 

9f . 
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i>j+2  _ 
9x9y 

9x 


where  i  =  1  through  N,  and  1  ^  j  £  M-2, 

Step  4 

The  equations  of  all  the  bicubic  polynomials  are  differentiated  to 

3  2  3  2 

give  9  s(x,y)/9x  9y  in  each  subrectangle.  Continuity  of  9  s/9x  9y  along 

each  horizontal  line  leads  to  a  set  of  simultaneous  equations  for  9f/9x9y 

at  all  the  interior  points  of  each  line.  The  necessary  values  9f/9x9y  at 

the  ends  of  each  horizontal  line  were  obtained  in  Step  3. 

The  equations  to  be  solved  can  be  written: 


where  1  £  1  £  N-2  and  1  £  j  £  M. 

Since  each  part  of  the  calculation  procedure  involves  the  solution  of 
simultaneous  tridiagonal  equations,  it  is  relatively  easy  to  compute  the 
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spline  by  Gaussian  elimination.  Moreover,  if  it  is  recognized  that  the 
tridiagonal  matrices  of  Steps  1  and  4,  as  well  as  those  of  Steps  2  and  3, 
are  identical,  the  number  of  Gauss  elimination  procedures  can  be  cut  in 
half.  (This  is  equivalent  to  the  statement  that  if  Ax  =  b  and  Ay  =  c  are 


to  be  solved  for  x  and  y,  then  A  need  be  determined  only  once.) 


Documen  tation 


This  section  documents  the  two-dimensional  spline  interpolating  pack¬ 
age  Implemented  in  Fortran  IV.  Subroutines  within  the  package  are  able  to 
return  interpolated  values  of  a  function  and  its  first  partial  derivatives. 
The  interpolating  functions  are  defined  within  a  domain  covered  by  a  rec¬ 
tangular  grid  in  the  x,y  plane.  Values  of  the  function  to  be  interpolated 
are  known  at  the  grid  intersections  on  and  within  the  boundary  of  the 
domain.  Values  of  the  first  derivatives  are  known  for  grij  intersections 


on  the  boundary,  and  values  of  the  second  derivative  3  f/3x^y  are  known  at 


the  four  corners  of  the  domain. 


Using  the  complete  interpolation  package  is  a  two-step  process. 

First,  the  user  creates  a  file  containing  raw  input  data  and  runs  the 
CREATE  program  to  create  a  binary  file  containing  spline  data.  Employing 
the  spline  data,  interpolating  functions  may  then  be  called  from  the  user's 
own  Fortran  programs  like  any  other  Fortran  real-valued  function.  The  raw 
input  data  file  must  contain  the  following  information  in  the  order  shown 
using  free- format ting  conventions; 


yi,y2.y3 


f(x^,y^),  fCx^.y^),  f(x^,y^) 


fCxj^.yf) 


fCxj.yj^)  ... 


3x  ••• 


3f  ,  . 


3x  •“ 


^ 

37 


\ 

Ty 


*-  *•*’  N*  ■*' 


Wi 


.*•  I.‘'  .  •  •  • 

n."  O  a'*  •  •  .  •*.  V 


SBBBKlfaJ 


••• 


^  (X  y  ) 

9y  ^  N’^M^ 


3x3y 


2  2 
'3x3y 


9^f  ,  . 


where  N  and  M  are  Integers  denoting  the  number  of  vertical  and  horizontal 
lines  in  the  grid.  Though  N  and  M  must  be  integers,  the  other  data  may 
take  any  legal  single-precision  real  value  (N  and  M  must  both  be  4). 

and  yj^  need  not  be  regularly  spaced;  the  grid  can  be  made  finer  in 
areas  where  greater  resolution  is  required. 

Two  additional  Inputs,  XYSC  and  FSCALE,  allow  for  data  scaling.  The 
scaling  parameters  allow  for  more  convenient  input  of  raw  data.  Their 
function  is  most  clearly  explained  by  Table  1.  If  XYSC  and  FSCALE  are 
omitted,  their  default  value  is  l.O. 

Execution  of  CREATE  generates  the  spline  data  file  subsequently  refer¬ 
enced  by  the  interpolation  functions.  The  CREATE  program  is  designed  to 
operate  through  the  intermediary  of  an  Interactive  terminal.  It  asks  for 
the  name  of  the  raw  data  file  and  the  name  of  an  output  file  in  which  to 
store  the  binary  spline  data  it  generates.  If  all  goes  well,  CREATE  will 
terminate  with  the  message  'normal  termination.’  Minor  code  revisions  can 
easily  be  made  to  alter  CREATE  input/output  functions.  The  CREATE  program 
listing  is  provided  in  the  appendix. 

The  interpolation  functions  SPLINE,  DSDX,  and  DSDY  return  the  value  of 
the  spline  S  and  its  first  partial  derivatives.  Calls  to  the  functions  may 
be  Incorporated  into  Fortran  code,  as  in  the  following  examples  evaluated 
a  t  x  =  XI ,  y  =  Y1 : 

HEIGHT  =  SPLINE  (XI  ,Y1  ,N) 

XSLOPE  =  DSDX  (X1,Y1,N) 

YSLOPE  =  DSDY  (XI ,Y1 ,N) 

GRADM  =  SORT  (DSDX  (X1,Y1,N)**2  +  DSDY  (XI ,Y1 ,N)**2) 


Each  of  the  Interpolation  functions  requires  the  spline  data  file  that 
CREATE  produced.  Befor'^  ^  Hing  any  of  the  functions,  spline  data  must  be 
read  into  a  COMMON  bio  single  call  to  subroutine  SPLOAD: 


CALL  SPLOAD  (<infi. 

where  <inflle>  is  the  Fortran  unit  number  assigned  to  the  spline  data  file. 
The  parameter  n  specifies  one  of  up  to  three  separate  spline  data  file 


I 


locations  in  COMMON,  l.e,  n  may  be  chosen  by  the  user  as  equal  to  either  1 , 
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Table  1.  Function  of  XYSC  and  FSCALE. 

Input  to  CREATE  =  factor  x  data  in  file 

X  =  XYSC  X  X 


y  =  XYSC  X  y 

^  _  FSCALE  ^  ^ 

3x  “  XYSC  9x 


^  ^  FSCALE  ^  ^ 

3y  XYSC  3y 


2,  or  3.  Thus,  the  same  program  code  can  call  upon  up  to  three  spline 
functions. 

SPLOAD,  SPLINE,  DSDX,  and  DSDY  are  listed  in  Appendix  A,  and  a  sample, 
called  PROGRAM,  demonstrates  the  use  of  these  routines. 


AVERAGE  LINEAR  VELOCITY  ROUTINE 


Description 


The  velocity  field  is  an  essential  input  to  the  calculation  of  fluid 
particle  trajectory.  For  so-called  laminar  flow  in  porous  media,  this 
information  is  supplied  by  the  average  linear  velocity,  that  is,  an  average 
velocity  of  the  particles  being  carried  along  through  the  pores: 


V*  =  -  JS  Vh 


where  K  is  the  medium  hydraulic  conductivity,  ifi  is  the  effective  porosity, 
and  Vh  is  the  gradient  of  hydraulic  head  h. 

Equation  20  presents  a  somewhat  idealized  view  of  fluid  flow  in  porous 
media.  The  actual  speed  of  particles  moving  through  the  void  spaces  ranges 
from  near  zero  at  the  grain  boundaries  to  many  times  the  magnitude  of  v*  in 
the  larger  spaces  between  grains.  This  phenomenon  accounts  for  some  of  the 
observed  smearing  of  sharp  solute  fronts  in  porous  media  flows.  It 
appears,  however,  that  the  effect  is  less  consequential  given  only  the 
approximate  descriptions  of  the  hydraulic  head,  hydraulic  conductivity,  and 
effective  porosity  fields  that  are  normally  available. 


Documentation 


Fortran  program  FG  is  written  to  supply  the  Cartesian  components  of  a 
two-dimensional  velocity  field  according  to  eq  20.  Used  in  conjunction 
with  a  program  called  MAIN,  FG  requires  spline  data  files  for  h  (n=l),  K 
(n=2),  and  (n=3).  Listings  of  programs  MAIN  and  FG  are  included  in 
Appendix  A. 


FLOW-LINE  DETERMINATION  ROUTINES 


Description 

The  equation  of  continuity  for  two-dimensional  flow  is: 


3x 


pv 


(21) 


where  p  is  fluid  density  and  t  is  time.  Using  the  characteristics  approach 
and  assuming  constant  density,  eq  21  becomes  the  linked  system  of  ordinary 
differential  equations: 


V 

X 


dx 

dt 


dp 

dt 


r  X 
=  -p[^ 


3v 

+  =  0. 

3y 


(22) 


The  first  two  equations  define  a  flow  line,  while  the  last  simply  states 
the  assumption  that  fluid  density  is  constant  along  a  flow  line.  The 
problem  of  determining  the  flow  line  becomes  one  of  solving  the  joint 
equations 


f (x,y ,t) 


dt 


g(x,y  ,t) 


for  x( t)  and  y( t) . 

Two  alternative  methods  for  solving  eq  23  are  proposed. 


(23) 


Me  thod  1 

A  fourth-order  Runge-Kutta  technique  is  used  to  start  the  procedure, 
which  is  continued  by  a  more  efficient  predictor-corrector  scheme.  The 
numerical  solution  method  starts  at  point  (xQ,yQ)  at  time  zero.  Using 
a  time  step  At,  successive  points  (xp,  yp)  along  the  trajectory  of  the 
particle,  which  began  at  point  (xq,  yo) ,  ate  obtained.  The  Runge-Kutta 
algorithm  for  accomplishing  this  is 
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(24) 


X  . ,  =  X  +  r  (s,  +  2a„  +  2a-  +  a.) 
n+1  n  6  1  2  3  4 

^n+l  ““  +  i  ^  2b2  +  2b3  +  b^) 

^  =  At  f(x^.y„,0 

a„  =  At  f(x  +  a,/2,  y  +  b  /2,  t  +  At/2) 
I  n  1  n  1  n 

b^  =  At  g(x^  +  a^/2,  y^  +  b^/2,  t^  +  At/2) 

a_  =  At  f(x  +  a-/2,  y  +  b  /2,  t  +  A  t/ 2) 

3  n  z  n  2  n 

b^  =  At  g(x^  +  a.^11,  y^  +  +  At/2) 

^  =  At  f(x^  +  a3.  y^  +  b3.  t^  +  At) 

b^=  At  g(x^+a3,  y„+b3,  t^+  At), 


(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 


and  f  and  g  are  defined  in  eq  23. 

Runge-Kutta  algorithms  belong  to  the  set  of  self-starting  numerical 
solution  methods.  Self-starting  means  that  the  determination  of  all  suc¬ 
cessive  points  (Xj^.yn)  requires  only  the  starting  point  (xQ,yQ).  In  other 
In  other  words,  calculation  of  Xj^  and  yj,  depends  only  on  the  known  values 
Xj^-i  and  yn-i«  The  set  of  nonself-starting  methods  requires  the  values 
(Xj,,yj^)  to  be  given  at  more  than  one  point  along  the  trajectory.  For  exam¬ 
ple,  a  fourth-order  predictor-corrector  algorithm,  called  Milne's  method, 
requires  the  values  xq,  x^,  X2,  X3,  yQ,  y^,  y2,  and  y3  to  calculate  suc¬ 
cessive  values  of  x^j  and  y^. 

In  addition  to  the  question  of  starting  values,  the  efficiency  of  the 
calculation  procedure  is  an  important  factor  in  selecting  a  numerical 
method.  It  turns  out  that  Milne's  algorithm  is  significantly  more  effi¬ 
cient  than  the  Runge-Kutta  method,  although  both  are  fourth-order  accurate. 
Method  1  takes  advantage  of  the  Runge-Kutta  self-starting  feature  and  the 
efficiency  of  Milne's  method.  Given  a  starting  point  (xQ,yQ) ,  the  Runge- 
Kutta  procedure  is  used  to  obtain  (xj,yj),  (x2,y2).  and  (x3,y3).  At  that 
stage  the  necessary  starting  values  are  available  for  Milne's  method,  which 
is  then  used  to  generate  succeeding  points. 
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Milne  s  predictor-corrector  method  consists  of  two  steps.  First, 
predicted  estimates  of  and  yn+i  are  calculated.  Let  these  be 

denoted  and  .  Second,  the  predicted  values  are  corrected  to 

obtain  the  final  values  and  y^+i  at  the  end  of  a  time  step.  The 

algorithm  is 


Then 


+  ^  [2x’  -  x’  +  2x' 
n-3  3  n  n-1  n-2 


n+1 


'’a-S  ^  ^  -  y;.,  .  2y;.^ 


n+1 


X  , ,  =  X 

n+1 


+  ^  fx*'  +  4x*  +  X*  1 

n-1  3  n+1  n  n-1 


■  'n-i  r  fCi  ‘‘K  *  K-i 


whe  re ; 


X*  -  f(x  ,y  ,t  ) 
n  n’^  n  n 


(34) 

(35) 

(36) 

(37) 

(38) 

(39) 


Consideration  of  the  Runge-Kutta  and  Milne  algorithms  shows  that 
Milne's  method  requires  only  two  evaluations  of  f  and  g  per  time  step, 
whereas  Runge-Kutta  requires  four.  This  makes  Milne's  method  more 
efficient. 

Nonself-s  tar  ting  methods  typically  assume  constant  At,  whereas  self- 
starting  methods  allow  for  change  of  At  at  each  time  step.  For  problems  in 
which  the  frequent  change  of  At  is  desirable,  exclusive  use  of  a  self¬ 
starting  method,  such  as  a  Runge-Kutta  algorithm,  is  advised. 

Method  2 

The  algorithms  of  Method  I  do  not  permit  the  stepwise  estimation  of 
accuracy,  nor  do  they  allow  for  time-step  modification  in  response  to 
algorithm  performance.  Method  2,  however,  incorporates  both  of  these 
fea  tures. 

Consider  the  following  Runge-Kutta  algorithm  from  Sarafyan  (1966); 

(4)  1 

+  r  ^3.  +  a  ) 


(40) 


=  X  +  ^  (l^a,  +  35a,  +  162a,  +  125a,)  (41) 

n+l  n  JJd  1456 

^ntl  =  +  3k  1^263  +  125b^)  (43) 

aj  =  At  f(x^.y„.t^)  (44) 

bj  =  At  gCx^.y^.O  (45) 

a^  =  At  f(x^  +  a^/2,  y^  +  b^/2,  t^  +  At/2)  (46) 

b^  =  At  g(x^  +  aj/2,  y^  +  b^/2,  t^  +  At/2)  (47) 

a^  =  At  f(x^  +  a^/4  +  a^/U,  y^  +  b^/4  +  h^/^,  t^  +  At/2)  (48) 

b^  =  At  g(x^  +  a^/4  +  a^/it,  y^  +  bj/4  +  b2/4  ,  t^  +  At/2)  (49) 

a^  =  At  f(x^  ~  ^2  ^  ^^3’  ~  ^2  ^  2b3.  t^  +  At)  (50) 

b^  =  At  g(x^  "  ^2  '*’  ^^3’  ~  ^2  2b3,  t^  +  At)  (51) 

a,  =  At  f(x  +  7a, /27  +  10a,,/27  +  a,/27, 

5  n  1  2  4 

y„  +  7b, /27  +  10b,,/27  +  b,/27,  t  +  2  At/3)  (52) 

n  1  2^0 

b3  =  At  g(x^  +  7a^/27  +  10a2/27  +  a^/27, 

y^  +  7bj/27  +  10b2/27  +  b^/27  ,  t^  +  2  At/3)  (53) 

f.  =  X  +  (28a,  -  125a„  +  546a^  +  54a,  -  378a,)  (54) 

n  625  12345 

n  =  y  +  (28b,  -  125b„  +  546b,,  +  54b,  -  378b,)  (55) 

n  625  12345 

a,  =  At  f(^,n,t  +  A  t/ 5)  (56) 

0  n 

b,  =  At  g(  0 ,1  ,t  +  A  t/ 5)  .  (5  7) 

b  n 


As  eqs  40  through  43  indicate,  the  Sarafyan  algorithm  is  an  embedded 
Runge-Kutta  scheme,  that  is,  the  same  information  used  to  calculate  the 
fifth-order  result  x^+j  can  be  used  to  calculate  an  embedded  fourth-order 

(4) 

result  with  essentially  no  extra  effort.  Wlien,  after  a  given  time 
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step  At,  x^+i  and  agree  to  m  significant  digits,  x^+i  nust  be  accurate 

to  that  number  of  digits.  This  feature  of  the  embedded  algorithm  is  used 
to  increase  or  decrease  the  time  step  so  as  to  balance  efficiency  and 
accuracy. 

Cartesian  coordinates  (x,y)  along  a  flow  line  are  calculated  at  points 
in  time  separated  by  time  steps.  The  degree  of  detail  desired  for  plotting 
a  flow  line  will  determine  a  selected  maximum  time  step  AT,  chosen  inde¬ 
pendently  of  accuracy  considerations.  For  each  interval  (t,  t  +  AT), 

Method  2  automatically  determines  a  constant  At  that  will  guarantee  that 
every  successive  result  within  the  interval  is  accurate  to  a  specified 
number  of  significant  digits  m.  The  rule  for  dividing  the  interval 
(t,  t  +  at)  requires  that 

AT  =  2^  At  (58) 

where  k  may  equal  0,  1,  2,  3,  or  4 .  When  accuracy  consistently  exceeds  nri-1 
significant  digits  throughout  a  time  interval  AT,  the  time  step  At  is  dou¬ 
bled  (but  never  exceeds  AT). 

Documentation 

The  approaches  of  both  Methods  1  and  2  are  implemented  in  the  form  of 
the  Fortran  subroutines  listed  in  Appendix  A.  The  two  subroutines  are  fully 
interchangeable,  and  operate  in  conjunction  with  the  MAIN  and  FG  programs. 

MAIN  serves  to  supply  either  method  with  necessary  inputs.  Execution 
of  MAIN  asks  the  user  to  give  the  file  names  of  spline  data  files  for  h,  K, 
and  <t),  as  previously  described.  It  then  asks  for  an  input  file  name  con¬ 
taining  data  on  time-step  size,  scale,  flow  region  boundaries,  and  origins 
for  the  flow  lines  to  be  calculated.  Finally,  MAIN  asks  for  a  file  name  in 
which  to  output  the  computed  points  along  the  flow  lines.  The  input  file 
has  the  following  free-format  structure: 


XYSC; 

xi ,  Xn,  yi ,  yn; 

’'0.  yo'> 


(or  AT  for  Method  2);  time  step 
scaling  factor  for  x  and  y  data 
flow  boundaries 
origin  of  flow  line  1 


’'0.  yo; 


origin  of  flow  line  Z 


•• 
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The  scale  factor  multiplies  all  XQ.yQ  data  and  flow  boundary  data  in 
the  file.  The  MAIN  program  is  designed  to  continue  reading  origins  until 
it  encounters  an  end-of-file.  Calls  from  MAIN  to  either  Method  1  or  2 
subroutines  include  passing  XQ.yg  and  At  (AT).  Either  subroutine 
outputs  flow-line  point  data  to  the  specified  file  at  intervals  At  (AT). 
Calls  to  subroutine  FG  originate  from  Method  1  or  2  routines. 

FLOW-LINE  PLOTTING 
Description 

A  plotting  routine,  listed  in  Appendix  A,  can  be  used  to  display  out¬ 
put  from  the  method  of  characteristics  solvers.  Plot  scale  is  set  by  the 
user  at  the  time  of  plot  execution.  This  scaling  feature,  used  with  trans¬ 
parent  drawing  films  or  tracing  papers,  enables  the  user  to  produce  flow¬ 
line  overlays  at  the  scale  of  an  appropriate  base  map  (see  Fig.  2). 

Documentation 

Program  SPLOT  is  designed  to  be  executed  interactively.  The  user  must 
provide  the  name  of  a  file  that  contains  output  from  either  method  of 
characteristics  and  a  plot  scale. 

Coordinate  data  stored  in  the  flow-line  file  are  given  at  tlie  actual 
base  map  scale,  i.e.  multiplied  by  scaling  factor  XYSC  (see  Documentation, 
above).  That  data  may  be  in  feet,  meters,  miles,  etc.  The  scale  factor 
for  plotting  (SCALE)  has  the  interpretation: 

1  inch  of  plotting  paper  =  (SCALE)  units  of  the  base  map 

If  coordinate  data  are  expressed  in  feet,  then  SCALE  =  2000  gives  a  plot  for 
which  1  in.  of  paper  equals  2000  f  t. 

Symbols  are  drawn  showing  the  location  of  a  fluid  particle  at  time 
intervals  whose  duration  is  printed  on  the  plot.  Tlie  user  can  choose  to 
have  symbols  drawn  at  every  Nth  point  of  the  flow  line  by  specifying  N;  all 
points  are  plotted,  but  only  every  Nth  is  given  a  symbol. 

APPLICATION 

The  flow  line  calculation  and  plotting  package  is  demonstrated  by  con¬ 
sidering  a  potential  groundwater  contamination  problem  at  the  hypothetical 
site  shown  in  Figure  5. 


Figure  5.  Base  map  for  demonstration  site. 


Figure  6  is  a  contour  map  of  steady  water-table  elevations  (in  feet 
above  MSL)  in  the  study  area.  A  water-table  contour  map  may  be  constructed 
by  interpreting  static  water  levels  at  a  number  of  wells,  by  observing 
topographic  and  stream  elevation  data,  and  by  considering  the  drainage 
behavior  of  the  study  region.  The  persistence  of  runoff  after  rainfall 
events  is  an  indicator  of  the  depth  to  the  water  table  in  the  vicinity  of 
streams;  sustained  baseflows  generally  indicate  a  high  water  table. 

Often  it  can  be  assumed  that  although  water  table  elevations  may 
gradually  change,  the  hydraulic  gradient  remains  fairly  constant  (steady 
state  or  quasi-steady  state).  This  may  be  a  valid  assumption  over 
considerable  lengths  of  time,  up  to  many  decades. 


Figure  6.  Water  table  contours  in  ft  above  MSL. 


For  the  purpose  of  demonstration,  hydraulic  conductivity  was  chosen  to 
vary  over  the  study  area  from  30  ft/day  in  the  south  to  50  ft/day  in  the 
north.  Effective  porosity  was  selected  as  0.1  throughout  the  study  area. 

In  actual  practice,  field  surveys  should  be  conducted  to  define  clearly  the 
spatial  variability  of  these  quantities. 

Raw  data  files  for  head,  conductivity,  and  effective  porosity  were 
prepared.  In  the  case  of  hydraulic  head,  data  were  specified  on  a  uniform 
grid  Ax  =  Ay  =  2000  ft.  Far  less  detailed  grids  were  required  to  describe 


K  and  5.  The  three  raw  data  files  were  analyzed  by  CREATE  to  generate 
spline  data  files. 


Figure  8.  Composite  flow  net  map. 


SUMMARY 

A  set  of  Fortran  computer  programs  for  the  calculation  of  flow  lines 
in  two-dimensional  velocity  fields  has  been  described,  documented,  and 
demonstrated.  The  approach  has  been  shown  to  have  useful  application  to 
the  prediction  of  contaminant  transport  in  groundwater  flow  systems. 
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APPENDIX  A:  PROGRAM  LISTINGS 
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*****  PROGRAI'^  CREATE  ***** 


This  routine  reads  a  file  of  raw  data*  analyses  it*  and 
generates  an  output  file  of  spHne  data. 


COMMON  /SPL/  F(5:.*Lj)» 
E'FDX(5C*5C1)» 
DFDY(^S0,50)* 
DFr}XY(5P,50) 
X  (  F  c  )  * 

Y  t 0  )  ♦ 

DELX (F  G  )  « 
DPLY(5'') 


/*  f<x»y) 

/ *  df/dx  (x»y) 

/*  df  /dy  (  X  *  y ) 

/*  d2f /oxdy  (x  fy ) 

/*  X  grid  values 
/*  y  arid  values 
/*  delta  X  along  grid 
/*  delta  y  alono  grid 


INTEGER*!  IN^ILE 
INTEGER*!  OUTFIL 


/»  logical  unit  #  of  inout  file 
/*  logical  unit  U  of  output  file 


FARAPETEP  infile  -  5  /«  only  file  unit  we  open 

PARAMETER  OUT'^IL  =  F. 


MNEMONIC  CODES  FOR  CRREL  FILE  CALLS 

INTEGER *5  JiREAD*dlWRlT,J$RDWR*J$CLOS*J$FRUD*JiRWND 


PARAMETER 

JIREAD 

- 

:i* 

/* 

READ 

JtUR  I  T 

Z, 

12* 

/  * 

WRIT  [ 

JSROwR 

z 

/* 

READ  AND  WRITE 

JICLOS 

z 

:  A* 

/* 

CLOSE  FILE 

JIDELE 

z 

:5* 

/★ 

DELETE  FILE 

J$FRW0 

z 

;f>* 

/* 

‘'OVE  FORWARD  ONE  FILE  (MAG  TAPE) 

J  J  R  <•  N 1' 

- 

:7 

/  * 

REWIND 

Open  input  fi 
Get  filename 

l  e 

from 

user 

call  GF  ills  ( 1  NF  ILt  »J<f.RE  AU* ’Sol  ine  data  infile»,lF.) 


&pt  data  fron  infile 


RLAutlNFlLE**) 
R  E AD (  INFI  LE  *  *  ) 
READ(INFILE,»> 
RE  AD ( INF  I LE  *  ♦  > 
READ ( INFILL  «*  ) 
«EAD( INFILL  ♦*  ) 
REaD(INFILE*») 
REAOdNFILE**) 
R  E  A  D  (  I  \  F  I  L  E  «  •  ) 


NX»NY  /*x*ydimerisions 

(x(I),l3ltNX)  /*  read  x  values 

(Y(I),I=1*NY)  /»  read  y  values 

( ( I  * J)  ♦  I =1 tNX  )  * J=  1  ,NY  )  /*  f(x*y)  at  Xnots 
(DFCX<1 « J)»J=1*NY)  /*  df/dx  at  left  side 

( DFDX ( NX , J) , Jrl tNY )  /*  df/dx  at  right  side 

(DFDYt I,  1),I=1,NX)  /*  df/dy  at  bottom 

(DEDY ( I *MY  )* I =1 *NX  )  /*  df/dy  at  top 

DFDXYtlfllfDFDXYtNX*!)*  /*  d?f/dxoy  at  corners 
DFDXY  < 1  *  NY) , DFDXYC  NX.NY ) 


Read  scalino  factors  from  the  file  (if  present) 
R  r.  A  0  <  I N  F  I  L  E  ,  ♦  *  E  N  n  =  A  )  X  Y  S  C 
GO  TO  1 
XYSC  =  1  .  0 

FrAD(INFILEf*»ENL=F)FSCALE 
GO  TO  5 
P  S  C  A  L  E  =  1  .  t 


Close  input  file 

CALL  FI  LE  SS  <  I  OF  IL  E  f  J  J  CL  0*'  ) 
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C  USE  COORDINATE  SCALE  FACTOR  TO  SCALE  INPUT  X»S  AN'  Y  ♦ ‘•' 

C  AND  FUNCTION  SCALE  FACTOR  FOR  F  AND  ITS  DERIVATIVES 

DO  2  1  =  1  ♦NX 

2  X(I )=X(I)*XYSC 
DO  3  I=1^NY 

3  Y <  I  > =Y( I ) *XYSr 

no  9  I=ltNX 
DO  9  J=lfNY 

9  F  (  I ♦ J )  =  F ( I ♦ J)  *FSC ALE 

DO  7  I=lfNX 

DFOY(Ifl)=DFCY(I,l)*FSCALE/XYSC 
7  OFOY(IfNY)  =  nFL'Y(IfNY)*FSCALE/XYSC 

DO  8  I  =  1  ♦  N  Y 

DFnXtl,I)=DFDX(l,I)*FSCALE/XYSC 
a  DFOX(NX,l>=CFnX(NX,I)*PSCALE/XYSC 

X YSC2  =  X  YSC*XYSC 

DFDXY<l^l)=DFDXY<ltl)*FSCALE/XYSC2 

DFDXY(1,NY)=DFDXY(1^NY)*FSCALE/XYSC2 

DFDXY(NX,1)=CFDXY(NX^1)*FSCALF/XYSC2 

;DFDXY(NXfNY)=nFDXY(NXfNY)*FSCALE/XYSC? 

C  Open  output  file 

C  Get  file  Hr  (re  frof\  user 

CALL  OF iLEi (0 JTF IL ♦ Ji«H 1 T ♦ ‘Out  put  file  name^tlO) 

C  Call  spline  to  generate  matrix  values 

CALL  SPLNXY(NX,NY) 

C  vrite  matrices  tc  output  tile 

•HITE  (OUTFIL)  NXfNY 

-RITE  (OUTfIL)  <  (F <  I  ♦  J  )  ♦  1  =!♦ NX  )  ♦ Jr 1,NY ) 

•RITE  (OUTFIL)  < ( DFO X < I ♦ J ) ♦ I = 1 ♦ NX > ♦ J= 1 ♦ N Y > 

WRITE  (OUTFIL)  ((OFOY(IfJ)^I=ltNX)fJ=l^NY) 

WRITE  (OUTFIL)  (  (D«^OXY  (  1  ♦J)  ♦  1  =  1,NX  )♦.)=!  ♦NY) 

WRITE  (OUTFIL)  ( X (  I  ) ♦  I  =  1  ♦  NX ) 

WRITE  (OUTFIL)  (Y(I)^I=ltNY) 

NX“1=NX-1 

NYM1=NY-1 

WRITF  (OUTFIL)  (DELX(l),Irl,V,Ml) 

WRITE  (OUTFIL)  (DELY(I),I=1^NYM1) 

C  Close  output  file 

CALL  FI LES$(CuTFILf JICLOS)  /♦  close  output  file 

C  Print  messaoc  to  ter'iiinal  indicating  norrnal  terrrination, 

CALL  TNOU  (•norsr.il  tertrtination’flR) 

999  CALL  EXIT 

TNO 


SUBROUTINF  3PlNXY(NX^  NY) 

C  PY  C.J.  DALY 

COMMON  /SPL/  rtSJ^bO)^  OFDX  (  b  :  ♦‘^  2  )  ♦  FIFO  Y(  SO  ♦  E  !?  )  ♦  FD  V  Y  (  S  C  ♦  5  0  )  ♦ 
1  X(bC)^  Y(bo),  UELX(^O)^  DFLY(S:. ) 


DIMENSION  TRIxU(bC)f  TRIXD(f):)^  TRIXL(S:) 


DIMENSIC'O  T«irb(bw),  T  K  I  Y  <  5  )  ♦  T-<!YL(^.C) 
0If*FNS10\  *5H?(50)«  DG.YXt'-:),  '31Y»‘'.-:) 


Compute  ccefficients  tor  the  snltne,  ':otr  that  Gausslar 
elimination  ir  cirrien  out  frr  the  triPianonal  matrices 
cetore  proceecin;  to  Steps  Ore  thrcu<)h  Four. 

\  <1^1  =  NX  -  1 

NY«1  3  NY  -  1 

N  >.  Y  2  r  NX  -  ; 

NY'i?  =  NY  -  2 

NXM3  =  NX  -  3 

NYf^.i  =  \1  -  }. 

Compute  cx'Lta  x  and  oelta  y . 

CO  1  1  =  1  ♦  N  V  M 1 

pfLX(I)  :  X(l4l)  -  X(I) 

CO  2  I  =  ItNYi-l 

HFLY  <  I  )  =  Y  (I  ♦!  )  -  Yd) 

Fill  tridiagonal  matrices  for  all  four  Steps. 

JO  3  I  =  1,NX''2 
TRIXUd)  =  CELXt  I) 

TSIXO(I)  =  <UtLXd>l)  ♦  DfLXtl))  *  2 

Tf-ixL  (I )  =  nrixt  1*1) 

DO  ^  I  =  l»NY-^-2 
TRl YU(I  )  =  DELY<  I  ) 

T=»IYD(I)  =  (DELYd*!)  •»  DCLYd))  *  2 
TRIYL(I)  =  CELY(dl) 

Perform  pivotina  of  triciaoonal  matrices^  downward. 

Save  pivot  values  for  all  four  Steps. 

PO  5  1  =  l,NX‘-'3 

TKIXL(I  +  1)  =  -TR  I  XL(  !♦)  )/TRI  XDd  ) 

ThIXD(I*1)  =  TRIXOd-*-!)  ♦  TRIXLCM1)*TRIXU(I) 

DO  b  I  =  ItNYx-S 

TRIYL(I  +  1)  =  -TP  I  YLd^l  ) /TRI  YD<  I  ) 

TRIYD(I-»1)  :  TKlYCd  +  l)  >  TR  I  Y  L  d  ♦  1  >  *  T  R  I  Y  U  (  I  ) 

Perform  pivoting  of  tridiagonal  matrices*  upward. 

Save  pivot  values  for  all  four  Steps. 
uO  7  I  =  1»NXY3 
d  =  NXM?  -  I 

TRIXUdI)  =  -TR  I  XU  ( I  I  > /TR  IXD  d  I +1  ) 

DO  a  I  =  1  1  N  Y  Y 3 

II  =  NYf'C  -  I 

TRIYUdl)  =  -TR  I  YU<  I  I  ) /TRI  YG  d  I  ♦!  ) 

Compute  ana  s*jve  ratios  ♦or  efficiency. 

CO  9  I  =  l*Ny"2 

DOTXd)  =  nELXd)/DLLX(l*l) 

L'D  i;  I  =  l*^Y^'^ 

POTYd)  :  DELYl  I  )/DrLY(  dl) 

Ferform  Step  One. 

DO  11  J  =  1  , N  Y 
DC  12  I  =  1*NXIX? 

RHS(I)  =  3.0  »  (  (  R  (  1  ♦  1  .  vJ> -F  (  I  ♦  J>  ) /DOTX  (  I  )  ♦ 

1  (F(  d2«J)-F(  I>1»J)  )*DOTX(  I>) 

CONTINUE 

PHSl’.'X.YD  =  RIS(CXM2)  -  1PIX"(NX*'2)*DFDX<NX»J) 

R  H  S  (  1  )  :  I-  H  S  d  )  -  T  P  I  X  L  d  )  *  D  F  [  ,  V  (1  »  J  ) 

0  C  13  I  =  1  *  N  X  P  3 
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13  RHS(I*1)  =  KHS<I*1)  ♦  TR1XL<  I+n*RHS(  I) 

DO  lA  I  =  ItNVMS 

II  =  NXM2  -  I 

lA  RHS(II)  =  RHS(II)  ♦  RHS < I !♦! ) *TR 1 XU( I  I  ) 

DO  15  I  =  2»\.<M1 

15  DFOXIUJ)  =  RHS  (I -1)/TR  I  XO(l -1  ) 

11  CONTINUE 

C  Perform  Step  Twc. 

DO  It  I  =  1«NX 

DO  17  J  =  1«NYM2 

KHS(J)  =  3.0  *  (  <F  (I  »  J+1 ) -F  (  I  tJ)  >/D(JT  Y  (J  »  ♦ 

1  (F(  I  ,  J  +  2)-F(  I.  >*OQTY(  J)  ) 

17  CONTINUF 

RHS<1>  =  RHSU)  -  TRI  YL  <1)*DFDY(1  *1) 

RHS<NYM2)  =  RHS(NYM2)  -  T R I Y U < NYP 2 ) *  OF D Y f I t N Y ) 
DO  IP  J  =  l,NYf13 

18  RHSIJ-*!)  =  RHD(J+1)  ♦  RHt  IJ)  ‘TkIYL  IJ*1  ) 

DO  19  J  =  1.NYM3 

JJ  =  NYM2  -  J 

19  RHS(JJ)  =  RHS<JJ)  ♦  PHS( JJ+1 )*TR1 YU(JJ> 

DO  2C  J  =  2tNYHl 

2C  DFDY(I.J)  =  RHS(J-l)  /  TRIYDU-D 
Ifc  CONTINUE 

C  Perform  Step  Three. 

DO  26  INDEX  =  1,2 
IF  (INDEX. EQ.l)  1=1 
IF  (INDEX. EO.D)  I=NX 
DO  21  J  =  1,NYM2 

21  PHS(U)  =  3.C  »  (  (OFDX(I,J*l)-r'FDXtl,J)  )/DOTY(J) 

1  (DFDXd  .J4.2)-DF0X(1  fJ^l))*DG:TY(J)) 

RHS(l)  =  RHS(l)  -  TRIYL (1 )*Of CXY(  I  ,1  ) 

R*-'S(NYM2)  =  RriS(.NYM2)  -  TR  1  Y  U  (  NYM  2  )  *DF  DX  Y  ( 1  ,  N  Y  ) 
DO  22  J  =  1,MYN3 

22  RHS(J^l)  =  RHS(J^l)  ♦  TRI YL(U+1)*RHS(J) 

GO  23  J  =  1,NYM3 

JJ  =  NYM2  -  J 

23  RKS(JJ)  =  RHSCJJ)  ♦  RHS( JJ+1 )*TRI YU( JJ> 

DO  2A  J  =  2,NYM1 

24  DFDXY(I,J)  =  RHS(J-1)/TRIYD(J-1> 

26  CONTINUE 

C  Perform  Step  Four. 

DO  2P  J  =  1,NY 

DC  2P  I  =  1,\>XN2 

29  FHSd)  =  3.0  *  (  (r)FDY(  I  ♦!  ,J> -LFDY  (  I,  J)  )/DGTX(  I  ) 

1  (DFDY(I*2,J)-DFDY(I>i,J))*ORTX(I)) 

RHS(l)  =  PHS(l)  -  TRTXLtl)*DFrXY(l,J> 

RHS(NXN2)  =  RHS(NXM2)  -  T R I  X U ( N XM 2  )  *DF D X Y ( N X , J > 
DC  30  I  =  1,NX,M3 

30  RHS(I+1)  =  RPS(I-H)  TRIXL(  I-*1)*RHS(I  ) 

DO  31  I  =  IfNXM? 

II  =  NX  M2  -  I 

31  RHS(II)  =  HHS(II)  ♦  TR I XU( 1 1 ) *RHS  (  I  !♦!  ) 

DO  32  I  =  2,NXM1 

32  nFDXY(I,J>  =  RHS ( I -1 ) /TR I XD ( 1 -1  ) 

28  CONTINUE 
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C  Spline  interpolating  function  package 


C  SPLOAD  *** 

SUPROUTINt  SPLOAD  (INFILEtK) 

C  Reads  spline  cata  into  one  of  three  (K  =  l»2t  or 

C  locations  in  COMMON. 

COMMON/SPL/  NX(3),NY<3)»  /•  number  of  x»y  indices 

1  r(fiO»5C*3>*  /*  function  values  at  interstices 

1  CFDX  ( 50 » ‘^0  *  3)  ♦  /»  dt/dx  at  interstices 

1  DFDY(50*50«5)*  /*  d^/dy  at  same 

1  nFDXY«50»5 0,3) »  /*  df/dxdy  at  same 

1  X(5C«3)t  /*  X  values  of  vertical  grid  lines 

1  YtSOtS)*  /»  y  values  of  hori?,  grid  lines 

1  lELX(50*3)*  /*  delta  x  alone  grid 

1  r.ELY(5C«3)  /*deltayalonn9rid 


INTEGER  INFILF  /*  data  file  for  spline  function 


REWIND  INFILE 

C  read  data  file  into  common 

1F{K,6T.3.0».K.LT,1)Kz1 

READ  (INFILE)  NX(K)»NY(K) 

KX  =  NX  (K  ) 

KY=NY<K> 

KXM1=KX-1 

KYMlrKY-1 

READ  (INFILE)  ( ( F ( I ♦ J tK ) ♦ i = 1 fK X > * J=1 tK Y ) 

READ  (INFILE)  ( ( DFDX ( I » d tK ) . I = 1 t K X ) * Jz 1 tK Y > 
READ  (INFILE)  ( ( DFD Y ( I ♦ J* K) ♦ I  =  1 • f X ) , J= 1  * K Y  ) 
READ  (INFILE)  ((DFDXY(I*JtK)»l=ltKX)*J=l.KY) 
READ  (INFILE)  (X(I*K),I=1«KX) 

READ  (INFILE)  (Y(I*K),I  =  1,'<Y> 

READ  (INFILL)  ( CEL  X ( 1 ,K  )  , I  =  1  ,  K XM 1  ) 

READ  (INFILE)  ( DEL Y ( I *K  )  , I = 1 »K YM 1  ) 

REWIND  INFILE 

999  REturm 
END 


C  ***  SPLINE  *** 


FUNCTION  SPLINE  (Xl,Yl,i')  /•  Returns  Real 

C  Function  spline(x»y)  returns  spline  interpolating 

C  function  value  at  (x»y). 

COMMON  /SFL/  NX(3)fNY(3)f  /«  numPer  rf  x.y  indices 

1  F(6CI»5Cf3)»  /*  function  values  at  interstices 

1  DFOX ( 50 t50 »3 ) «  /*  nf/dx  at  interstices 

1  uFDY(50«50f3)*  /*df/dydtsame 

1  ;FDXYI5 : »5C»3) ♦  /*  df/dxdy  at  same 

1  /•(5C;*3)»  /•  X  values  of  vertical  grid  lines 


SPLINE 


Paq#* 


1 

Y  (  5  D  »  3  1  « 

1 

DELX (5C»3)  * 

1 

CEL Y( 50.3) 

RFAL  X1,Y1 

/  • 

INTEGER  I»J 

♦  y  values  cf  horiz.  ria  lines 

•  Celta  X  along  qrid 

*  relta  y  along  qrid 

i nt erpo la t inq  tuncticr  arpuments 
su" r ec t ano I e  inoices 


C  Define  interpolating  furctions 

CX(X1*I*K)  =  (  ( ( X  1-X t I ♦! ,K ) ) **2) /DEL  X n tK ) * *? > ♦ 

1  <?*(^Xl-X(I*^))*<Xi-X<l■»l♦'))**2)/(DELX(I,K)**?)) 

cxi(xi*i»K)  =  (((xi-x(ut'’))**r>/[irLxii»K)**?)- 
1  (2»((Xl-X(I*l,K))*tXl-X{I»K))**2)/(DELX(I,K>*^M) 

CY(Y1»I«K)  =  <((Y1-Y(M1»*<)>*»2)/l'FLY<I»K)**2)-» 

1  (2»<(Yl-Y(ItK))*<Y]-y(I*l,X))**2>/(rELY(l,K)**3)) 

rYl(Yl,ItK)  =  <((Yl-Y(I,Kn**2)/DrLY(I«K)**2)- 
1  (2»((Y1-Y<I^1»K))*{Y1-Y(I,><))»*2)/(DELY(I,K)**1>) 

DX«X1,I»K)  =  (Xl-X ( 1 «K )  )  *  (Xl-X C 1*1 ,K ) ) **2/DELX (1 tK )* ♦? 

DX1(X1*I*K)  =  «X1-X(I^1,K))*(X1-X(U‘<))**2/0ELX(1,K>**2 

DY(Y1«I»K)  =  (Yl-Y(I»K))*tYl-YII+l,K))**2/DELY(l«K)**? 

DYKYltlfK)  =  (Yl-Y<I+ltKl)*<Yl-Y(ltK))**2/DFLYtI»K)**2 

C  *'ake  sure  the  point  is  on  the  map 

rF(K.GT.3,0R.K.LT.l)K=l 

KXrNX (K ) 

K Y  =  NY (K  ) 

IF  ( Xl.LT  .X <1  ♦K ) .OR.Xl .01 .X(KX,K ) )  GOTO  777 

IF  ( Y1  .LT  .Y  (1  ♦K ) .OR .Y1 .GT .YiKY.K ) )  GOTO  777 

C  Determine  I  and  J 

NXMl  =  KX  -  1 

NYMl  =  KY  -  1 

DO  13  I  =  I.NaRI 

IF  (xi.le.x(i  +  i.k))  go  to  if, 

C  exit  with  Xd.K)  <=  XI  <=  X(I  +  1,K) 

IQ  CONTINUE 

15  no  2.3  J  =  1»NYM1 

IF  (Yl.LF.Y(J*l.K))  GOTO  33 

C  exit  with  Y(J.K)  <=  Y1  <=  YiJ'el.K) 

20  CONTINUE 

C  Calculate  spline  interpolate 

3D  SPLINE  =  F(I.J«K)*CX<X1,I.K)*CY(Y1,J.>^)  ♦ 

1  F ( I , J+1 .K ) *CX< XI  .  I  ,K)  *C Y1  { Y1  ♦ J.K)  ♦ 

1  F(T  +  1,J,K)*CX,  l(Xl,l,K)*rY(Yl»J«K)  ♦ 

1  F(I>1,.J+1.K)*CX1«X1.I.K)*CY1(Y]»J.K)  ♦ 

1  nFOX(l.J»K)*DX<Xl,T,K)*CYfYl.J,K)  ♦ 

1  DFDX ( I  . J4l »K ) ♦CX  (X 1 ♦ I  .K ) *0 Y1  < Y1  . J »K  )  ♦ 

1  OF  OX ( I *1 , J.K) *  DX  1  <  XI  .  I  ,  K) *C Y( Y1 , J. K)  ♦ 

1  UFnX(I*l,j4l.K.  )*rxi<Xl,I,K)*CYl(Yl»J.K)  ♦ 

1  DFDY ( 1  ♦  J.K  )  *CX < X  1  .  I  ,K  )  *DY ( Yl. J ,K )  ♦ 
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I  DfDY ( I «J41 tK )»CX CXI,  I  .K)  *DY1 <yi , J,K)  ♦ 

1  DFDY(I*1,J,K>*LXHX1,1  ,K)*DYCY1,J,K)  ♦ 

1  :^FDY(l+l,j4l,K)*C>l(Xl,I,K)*ny!  (Y1  ,J,K)  ♦ 

1  DFDXYC I,  J,K) *0X( X1,I, K) *DY ( Yl, J,K)  ♦ 

1  DFDXY(1,J41*K)*0X(X1,1,K)*0V1(Y1,J,K)  ♦ 

1  DFDXY(l4l,J,K)*DXHXl,I*K)*0YCY!,J,K)  * 

1  DFDXY(!4l,j4l,K)*rxl(Xl,I,KI*DYl(Yl,J,K) 


999  RFTUKN 


C  error:  here  If  point  out  of  hounds 

777  WRITE  (1,1013)  XI, Y1 

1010  FORMA  T(  ‘.SPL  INF  ARGUMENTS  OUT  OF  ECUNOS:  •  ,  1  P  2  F  1  0  .  fe  ) 
STOP 


END 


C  ***  OS/DX  *** 


FUNCTION  DSPX  {X1,Y1,*') 


/*  returns  real 


C  Returns  the  partial  w.r.t 


X  of  the  spline  interpolating  function 


COMMON 

/SPL/  NX(3),NY(3), 

/* 

1 

F(5C,S0,3)  , 

/* 

1 

nFDX{bO,bfl,3), 

/* 

1 

CFOY(50,hO,3), 

f* 

1 

DFDXY<5G,5C',3)  , 

/* 

1 

X(50,3) , 

/* 

1 

Y(50,3), 

/♦ 

1 

DELX(50,3) , 

/* 

1 

DELY(50f 3) 

/♦ 

number  of  x,y  indices 
function  values  at  interstices 
df/dx  at  interstices 
Jf/dy  at  same 
df/dxdy  at  same 

X  values  of  vertical  grid  lines 
y  values  of  horiz.  grid  lines 
aelta  x  alonn  grid 
delta  y  along  grid 


REAL  X 1  ,  Y  1 
INTEGER  I,J 


/•  function  arguments 
/*  subrectangle  indices 


C  Define  interpclat ing  sub-functions  (  o  for  partials) 

PCX<X1,I,K)  2*(Xl-XCl4l,K)  )/nELXCI,K)**2  ♦ 

1  2*(X1-X(I41,K))**?/DELX(I,K)**3  ♦ 

1  9*{X1-X<I,K))*(X1-X(I41,K))/0ELX(I,K)**3 

PCXlfXl,I,K)  =  2*(X1-X(1,K))/DELX(I,K)**2  - 
1  24CXl-X(I,K))**2/nrLX(  T,K)**3  - 

1  R*<xl-XCI,K) )*(X1 -X (I  41 ,K ) )/DELX (1 ,K ) 4*3 

PDX(X1,I,K)  =  2*(Xl-X(I,K'))*<Xl-X(l4l,K))/DELX(I,K)**?  4 
1  (Xl-X( T4l,K)>  442/DELX(  1  ,K)*42 

PUX1<X1,I,K)  =  2* CX1-X( I  ,K) ) * ( Xl-XCI 41 ,K > >/DELXt I ,K) 4*2  4 
1  ( Xl-X (I ,K) ) 4*2/DELy ( I  ♦K) *4? 


CY(Y1,I,K)  =  ( ( ( Yl-Y< l4l,K) ) **r)/OELYC I,K)**2) 4 

1  (2*<(Yl-Y(r,K))4(Yl-Y(l4l,K))**2)/(DrLY(I,K)**3)) 
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CY1(Y1*I*K)  =  (<  (Yl-Y  ( I  «*<>)**?  >/DTLY(  I  tK  )  **2>- 
1  (2•((Y1-Y(I♦1♦K))*(Y1-Y(1♦K)>**2>/(DELY(I♦K)**:^)) 

DY(Y1,I,K)  =  tYl-Y<I*K)  )«(Yl-YCI^l,K)>**2/DTLYUfK)**2 
DY1(Y1«I»K)  =  <  Yl-Y  (T  ■»!  ♦K  )  )*  <  Yl-Y  (  1  *K  )  )**2/DEL  Y(  I  fK  )  **2 

C  Wake  sure  the  point  is  on  the  B.ap 

IF(K.GT.3,0R.h.LT.l)K=l 

KX=NX  (K  ) 

K YzNY  (K  ) 

IF  (X1.LT.X<1*K).0R.X1.G7.XIKX,K))  GOTO  777 
IF  (Y1.LT.Y(1*K).0R.Y1.GT.Y<KY,K))  GOTO  777 

C  Determine  I  aro  J 

NXMl  =  KX  -  1 

MYPl  =  KY  -  1 

DO  1C  I  =  1*NXM1 

IF  (X1.LE.X( 1+1»K))  GOTO  15 

C  exit  with  X<I»K)  <=  XI  <=  X(1>1,K) 

10  CONTINUE 

15  DO  2C  JzliNYMl 

IF  (Yl.LE.  Y<  J-i'l.K))  GOTO  3C 

C  exit  with  Y(J*K)  <=  Y1  <=  YCJ*1*K) 

20  CONTINUE 

C  Calculate  ds/dx  interpolate 

30  DSDX  =  F(ltJ*K)*PCXlXltl.K)*CY{Yl*J,K>  ♦ 

1  Fri,J*l»K)*PCX{Xl»I*K)*CYl(Yl»J,K)  ♦ 

1  F(I*1,J«K)*PCX1CX1*I*K)*CYCY1*J*K)  ♦ 

1  F(l4l,J>l,K)*PCXl(XltlfK>*CYl(Yl»JtK)  ♦ 

1  DFDXn*JtK)*PDX(Xl*ItK)*CYlYltJ*K)  ♦ 

1  DFDXd  tJ*l.K)*POX(Xl»  I,K)*CY1  (Y1»J*K)  ♦ 

1  DFDX <I JfK) *PDX1 CXI *I tK) *CY tYl* J»K )  ♦ 

1  DFOX(I»l*J*lfK)*PDXlCXlfI»K)*CYl(Yl,J*K)  ♦ 

1  DFOY ( I ,dtK) «PCX(X1,I,K) *OYI Yl, J*K>  ♦ 

1  DFOYCI »J*1*K)*PCX(X1*I*K)*DY1(Y1*J»K)  ♦ 

1  DFDYCI*l,J,K)*PCXl(Xl,l,K)*DY(Yl*d*K)  ♦ 

1  DFDY(I-fl»J  +  ltK)*PCXlCXl*I*K)»DYl(Yl,JtK)  ♦ 

1  DFDXY(],J,K)*PnX(Xl*I,K)*DY(Yl,  JfK) 

1  OFDXY(I,J*l*K)*PnX(Xl,I,K)*DYl(Yl»J*K)  ♦ 

1  DFDXYCI-»l.J*K)*PDXl(Xl,l*K)*DY(Yl»JfK)  ♦ 

1  DFDXY(I*l.J^l,K)*PDXlCXltI*K)*DYl(YlfJ»K) 

999  RETURN 

C  error;  here  it  point  out  oT  bounds 

777  WRITE  (IflClC)  X1»Y1 

ICIO  FORMATC’DSDX  ARGUMENTS  OUT  OF  BOUNDS;  •♦1P2E1G,6) 
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FUNCTION  OSDY  (XliYl*K) 


t*  returns  real 


‘Upturns  tte  partial  w.r.t.  Y  of  the  spline  Interpolating  function 


COMMON  /SPL/ 


NX(3) ♦NYt? ) , 

F  ( 0  ♦  5  0  *  3  )  . 
CFOXlSOtSO,?), 
nFDY(50*5C  t3)* 
0FDXY{5Ct5Q»3)» 
X(5C*3) » 

Y  (  5  :  t  3  )  ♦ 
:jFL''(5C*3)  • 
P£LY(50»3) 


number  of  xty  Indices 
function  values  at  Interstices 
df/dx  at  Interstices 
df/dy  at  same 
df/dxdy  at  same 

X  values  of  vertical  grid  lines 
y  values  of  horlz.  arid  lines 
delta  X  alone  grid 
eelta  y  along  grlo 


REAL  X1«Y1 


/»  function  arguments 


INTEGER  I»J 


/*  subrectannle  Indices 


Define  Interpolating  sub-functions  (p  for  partlals) 


PCYIYI.UK)  =  2*  <  Yl-Y(  I  ♦IfK)  )/DELY(I,K>**2  ♦ 

1  2* f Yl-Y( I^1«K) ) **2/DEL Y( I,K)**3  ♦ 

1  4*(Y1-Y11*K))*(Y1-Yn  +  1»K  )  T/DEL  Y  11  ,K  )  **3 


PCY1<Y1»1*K)  =  2*  <  Yl-Y  <  I  «K)  ) /;JELYI  I,K)**2  - 
1  2*(Y1-Y(I,K)>**2/DELY< I»K)**3  - 

1  4*(Y1-Y(I,K)  )*(Yl-Y(I-el,K))/DELY(I*K)**3 


PnY(Yl*I*K)  =  2*{V1-Y(I fK))*(Yl-Y<I*l*K))/0ELY(IfK)**2  ♦ 

1  ( Y1 -Y ( I ♦! ♦K ) )**2/DLLY( I fK)**2 

PCY1<Y1*I«K)  =  2* ( Yl-Y( I,K) ) *{ Yl-YC !♦! ,K) )/DELY( I ,K) **2  ♦ 
1  ( Yl-Y (1 tK) )**2/DLLY(I .K )**2 


r.X(Xl,I,K)  =  (  (  (  Xl-X(  !♦  1  »K)  )  **2I/0ELX(  I  »K)**2)  ♦ 

1  (2»(<Xl-X<T*K>)*(Xl-X(I^l,K))**2>/fDELXfI,K)**3>> 


CXKXl.IfK)  =  (  (  (  Xl-X  (  I  »K  )  >**2  )/DELX(  I  »K  )  **2)- 
1  (2*((X1-X(I+1»K))*«X1-XCI»K))*«2)/(DELX{I»K)**3)) 


DX(Xl,IfK)  =  (X1-X(I «K )  )*(X1-X ) )**2/0ELX(l ,K )* *2 

DXKXl.I.K)  =  (Xl-X(l -H  fK  >)♦  (Xl-X  CI,K  )  )**2/DELX  (  I  fK)  **2 


'^ake  sure  the  point  1s  on  the  map 
IF(K.GT.5.0R.K,LT.1)K-1 


X=NX  (K  } 

Y  =  NY  (K  ) 

IF  < XI .LT .X ( 1 ,K) .OK .XI , GT .X< KX,K) >  GOTO  777 
IF  (Y1.LT.Y(1»K).0R.Y1.CT.Y(KY»K>)  GOTO  777 


Determine  I  a-td  j 


NXMl  =  KX  -  1 
NYMl  =  KY  -  1 


DO  1 :  I  =  1  ♦  w Y M 1 

ir  (Xl.LF  .X  (  I-fl  ♦K  )  )  GOTO  ID 
exit  with  XU»K)  <=  a1  <=  X(1+1,K) 
rONTlNUE 


Itt:-  ■,;< 


v'*.- 


•.'.•■.•■■a 

***  -  Xj  _  , 
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15  DO  20  J=1.NYK1 

IF  ( Y1 .LE. Y ( J+1 *K  )  )  GOTO  3" 

C  exit  with  Y(dtK)  <=  Y1  <:  Y(J+1,K) 

2C  CONTINUE 

C  Calculate  ds/oy  interpolate 

30  DSDY  =  F(  It J«K)*f XlXl.I *K)*PCY(Y1*J*K)  ♦ 

1  F(I,U-*-l,K)*CX<XltItK)»PCYllYltJtK)  ♦ 

1  F(I*lt  JtK)*CXl<XltItK  )*PCY{YltUfK) 

1  F  (I  ♦!  t  J-»l  tK  )*CX1  iXltl  tK  )  *PCY1  (Y1  td  tK) 

1  DFDX(I,d,K)*DX(XltI,K)*PCY(Yl,dtK)  ♦ 

1  DFDX(I,d-»ltK)*DX(XltltK>*PCYl(YltdtK)  ♦ 

1  OFDX<I*ltdtK)*DXl(Xlt ItK)*PCY(YltdtK)  ♦ 

1  DFDX(I*ltd*ltK)*DXl<XltltK)*PCYl(YltdtK)  ♦ 

1  DFOYd  td,K)*CX<XltItK>*PDY(Yl,dtK)  ♦ 

1  DFDY  1 1  td-f]  tK)  ♦CX  <X1,  I  ,K)  *PnYl  <  YltdtK)  ♦ 

1  DFOY  ( I -t-l  tdtK  )  *CX1  (Xlt  I  tK)  *PnY(  YltdtK  )  ♦ 

1  DFOY(I-»ltd*ltK)*CXl(XltItK)*PDYlfYltdtK)  ♦ 

1  OFDX  Y<  I  tdtK)  *DX<  XI 1 1  tK)  *PDY  (  Y1 1  dtK  )  ■» 

1  OFC'XYt!td-»ltK)*DX<XltI»K)*FnYl(YltdtK>  ♦ 

1  OFnx Y ( 1+1 tdtK) *0X1 (XI t 1 tK) *PDY ( Y1 tdtK)  ♦ 

1  DFUXY(T-»ltd*ltK)*DXl(XltItK)*PDYl(YltdtK) 

999  RETURN 

C  error;  here  if  point  out  of  bounds 

777  WRITE  (1,1010)  XltYl 

ICIO  FORMAT(*nsOY  ARGUMENTS  OUT  OF  BOUNDS:  •tlP2E16.£.) 
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*****  Pt<CGRAP‘  ***** 

FiampLe  program  for  Spline  I  n  t  c  rpo  I  at  i  r»g  Package 

Tt  returns  function  value  anc  gradient  at  points  chosen  by  user 
MNFMONIC  CGPFS  FOR  CR"tL  FILF  CALLS 

INTEGFK  *?  JSRFAD*  JSyPIT  »  JiRDlrP  ♦JICLOS*  J$F-  LIDtJSRWND 


P ARAWFTEP 


1 

J iREAu 

z 

lit 

/* 

1 

JiWRIT 

z 

:2t 

/* 

1 

JSRDWR 

z 

:3, 

/* 

i 

JSCLQS 

z 

:^t 

/* 

1 

J JDELE 

z 

:Gt 

/* 

1 

JSFR JO 

z 

:6t 

/* 

1 

JiRWND 

z 

:7 

/* 

READ 

WRITE 

READ  AND  WRITI 
CLOSE  PILE 
DELETF  FILE 
MOVE  FORWARD  ONE  FILE 
PEW  1  NO 


(MAG  TAPE) 


INTEGER  INFILf 


/*  fortran  unit  number  of  data  file 


REAL 

Xt  Y* 

/* 

1 

HEIGHT  * 

/  * 

1 

XSLOPEt 

/* 

i 

YSLOPE 

/  * 

print  we  pass  to  spl 1ne 
interpolated  function  value 
ds/dx  at  X  *  y 
ds/dy  at  x  *y 


open  the  data  file  maoc  by  CREATE 

IMFILE  =  5  /*ft>rtranunit  JJofdataflle 


Get  ■filename  from  user. 

CALL  GFILEJ  (I  ME  ILE  ♦JiRE  AD»*SPLIfvlL  FILE**11) 


load  the  spline  data  file 

Data  is  arbitrarily  loaded  Into  location  K=1  in  COMMON. 
CALL  SPL0AD(I\FILE*1) 


close  data  file 

CALL  FILESt.(INFILEtJ$CLOS) 

get  point  frofr  user 

CALL  TNO'JA(*point  to  compute 

READ  (1,*)  X,Y 


call  interpolating  functions 


HEIGHT 

=  SPLINE ( X tYt  1  ) 

/» 

find 

function  value  at 

XSLOPE 

=  DSnXfXt Ytl) 

/♦ 

f  ino 

as/dx 

YSLOPE 

=  DSDY(XtYtl) 

/* 

find 

ds/dy 

PRINT  Ij)  HEI GMT »XSL0PE  .YSLOPt 
F0RMAT(1X,F1':'.?,2F1C.‘1) 

print  answers  on  terminal 


GOTO  10  /*  loop  until  user  gets  bored 
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C  Main  pronram  for  2-0  particle  movement  In  a  flow  field. 

C 

C  MMFMOMC  COOP';  FOR  CRREL  FILE  CALLS 

IME&ER*2  JlREADi  JSURIT  ♦JSROwR  »  J$  CLOS  *  JS  F  w  U  D  .  J  tR  UM  C 


PARAMETER 


1 

JSREAD 

z 

:  1  * 

/* 

1 

JS'JRIT 

r 

:2* 

/* 

1 

JSR  OWR 

r 

:3* 

/* 

1 

JSCLOS 

z 

:4  ♦ 

/* 

1 

JSDELE 

z 

:5* 

/* 

1 

JSFRwD 

z 

:6* 

/* 

1 

JSR  UNO 

z 

:7 

/* 

READ 

WRITE 

READ  AND  wRITF 

CLOSE  FILE 

DELETE  FILE 

MOVE  FORWARD  ONE  FILE 

REWIND 


(MAG  T  APE ) 


INTEGER  SPFILL»OUTPUT 

C  Points  along  a  single  flow  line  (XP*YP) 

niMENsiON  xp(  icoa)«YP(ii:oO> 

C  Boundaries  of  the  flow  region 

COMMON/ BOUND/ XI ♦XN«Y1 »YM 
I NF ILE=5 
OUTPUTrG 
SPFILE=7 
PRINT  lo: 

lOu  F0RMAT(//1X»4HG1VE) 

C  Read  hydraulic  head  spline  data  file. 

CALL  GFILE$(SPFlLE*J$READt*H  SPLINE  INFILE*»15) 

CALL  SPLOADlSPFlLEtl) 

CALL  FI LES$(SFFILE» JICLOS) 

PRINT  100 

C  Read  hydraulic  conductivity  spline  data  file. 

CALL  GFILE$(SPF  ILE»JwREAD,*K  SPLINE  INFlLE*»lEi) 

CALL  SPL0AD(SPFILEt2) 

CALL  FI LES$(SHF1LE«J$CL0S) 

PRINT  i;c 

C  Read  effective  porosity  spline  data  file. 

call  GFILESlSt’FlLE  tJSPEADt»PSI  SPLINE  INF1LE»*17) 

CALL  SPL0AD<SPFILE.3) 

CALL  FI  LESi(SFF ILEt JSCLOSI 
PRINT  100 

C  r-et  filename  for  file  containing  delta  t*  scale*  flow 

C  region  boundaries*  and  initial  points  of  flow  linrs. 

CALL  GFILE1 (I^FILE*J1KrAC♦•POI^T  INPUT  FILF»*1G) 

PRINT  ICO 

C  Get  filename  for  output  of  flow  line  points. 

CALL  GFILEKOuTPUT  *JWwRlTv«STPEtMLlNt  OUTPUT  FILE»*?2) 

PEAD( INFILE** )DELT 
READl INFILE ♦*  )XYSC 
READ(INFILE**)X1*XN*Y1*YM 

C  Scale  Input  flow  boundary  data. 

X1=X1 ‘XYSC 
X  N=  X  N  *  X  Y  S  C 
Y1=Y1 *XY3C 
YM=YM*XYSC 
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WRITE  (OUTPUT, 1000) OELT 
WRITE  (OUT PUT,lC3()Xl,XN,Yl,yp 
1000  FORMAT( 1X,AF1 0.2) 

Read  (X»Y),  the  beginninti  point  of  a  floy  line. 
Scale  input  values. 

READ(  INFILE,*,ENU=1)X,Y 
X=X*X YSC 
Y=Y*X YSC 

C  Compute  floy  line  coordinates. 

CALL  STREAM(OUTPUT,X,Y,r;ELT) 

C  Co  back  to  reaa  initial  point  for  another  flow  line 

GO  TO  2 

CALL  FI  LESidNFILEfJiCLOS) 

CALL  FI  LESS (OUTPUT, JSCLOS) 
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C  Average  linear  velocity  calculation  routine. 

SUr^ROUTINE  FG(T»X,Y»rfG,IOUT  ) 

REAL  K 

K  =  h,  ydraulic  conductivity  (isotropic).  C.O"'^*ON  location 
H  =  hydraulic  head.  C-pline  data  in  COMMON  location  1. 
PSI  -  effective  porosity.  COMMON  location  3. 

The  flow  region  is  XI  <  X  <  and  Y1  <  Y  <  Ym, 

COMMON/BGUND/Xl  ♦XNfYl  ♦YM 

Check  to  see  if  >  ana  Y  are  withirt  the  flew  recion. 
1F(X.LT.X1.0R.X.GT.XN)GO  TO  1 
IF(Y.LT.Y1.0R.Y.GT.YM)G0  TO  1 

Use  spline  files  to  compute  the  following. 

(DHDXtDHDY)  is  the  hydraulic  gradient. 
nHDX=DSDX (XtY«l) 

DhDY=DSDY(X« Y,l) 

K=SPLINE(Xt Y»?) 

PSI=SPLINE(X»Y*3) 

Darcy’s  law  for  average  linear  velocities. 

F  =  X  component  of  the  average  linear  velocity. 

G  =  y  comoonent  of  the  average  linear  velocity. 

F  =  -(K  *DHDX) /PSI 
G=-(KsOHnY) /PSI 
RFTUR  t: 

C  lOUT  =  1  sionals  the  crossing  of  a  flow  boundary. 

1  lOUTrl 

R  1 T  U  R  rj 


oonon 
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*****  METHOD  1  ALGORITHM  ***** 
SUBROUTINE  STREA''(OUTPUTtX«Y*()ELT) 


SOLVES  F(X»YtT)  CC/DX  ♦  6(X,Y«T)  DC/DY  ♦  DC/DT  =  0 

BY  THE  METHOD  OF  C  H  A  K  AC  T  E  P  I  OT  1  C  S  USl^iG  FOURTH  ORDER 
RUNGE-KUTT A  AND  HR E D I C T OR - CORR E CT 0 R  ALGORITHMS. 


INTEGER  OUTPUT 

DIMENSION  XNPM ( 3 ) , YNPM( 5 ) ,XNT( 1 OCl ) . YNT ( 1 COl ) 

IOUT=: 

C  when  IOUT  is  kETURNED  AS  1*  THE  CHARACTERISTIC  WILL  'E  CROSSING 

C  A  BOUNDARY  OF  THE  FLOW  REGION. 

NPOINT=l 
DELT2=DELT /2. 

D£LT^=DELT/?. 

DELTG=0ELT/6. 

DELT^3=A.*DEL  T3 
DELT83=8. *OELT3 

C  INITIALIZE  AT  THE  POINT  OF  ORIGIN  OF  THE  CHARACTERISTIC 

X  N  T  (  1  )  =  X 
YNT  < 1  Y 
XN  =  X 
YN  =  Y 
T\=0.0 

CALL  FG(TN*XN»YN*F1*G1*  IOUT) 

IF (IOUT. NC. 2)00  TO  1 

C  PERFORM  RUNGE-KUTTA  FOR  THE  FIRST  THREE  POINTS  TO  GF T  STARTED 

C  CLASSICAL  (4**>  FORMULATION. 

00  2  IT=1.3 
TN1=TN+C£LT 
XVAL  =  XN  +  Fli*DELT2 
YVAL  =  YN-»G1*CElT2 
TVAL  =  TN-*-0ELT2 

CALL  FG(TVAL*XVAL«YVAL»F2»G2*I0llT) 

IFdOUT.NE.OGC  TO  1 

XVAL  =  XN-»F2*DELT2 
YVAL=YN+u<'*DLLT2 

CALL  FG(TVAL,XVAL,YVAL»F3»G3.ir.UT> 

IF ( IOUT .NE. C) 00  TO  1 

XVAL=XN+F3*DELT 
YV  AL  =  YN-»G5*DELT 

CALL  FG  (TNI  .X  V  AL  *  YVAL  f  F  GA  .  i  OUT  ) 

I F ( IOUT  .UE  .  n ) oO  TO  1 

XNl  =  XN*DLLTG*{Fl->FA)-*r}ELT3»(R2-»F3) 

YM  =  YN*DELTF»(G1*GA)+LELT3*(C2+G3) 

XNT  (  I  T+  1 )  =XM 
YNT ( I T*1 ) =YN1 
NF0INT=IT+1 
XN=XK1 
YN  =  YN  1 
T  N  =  T  N  1 

CALL  FG(TN.XN.YNfFl*Gl.ir:uT) 

IFdOUT.NE.OGO  TO  1 

XNPM( IT)=Fi 
2  YNPM(IT)-01 


C  USE  FREDl CTOR -CORKFCTOR  METHOD  FOR  THE  REST  OF  THF  CHARACTERISTIC 

C  STANDARD  FOURTH  ORDER  MILNE  FORMULATION. 
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DO  3  IT=4,1jC1. 

TM  =  TN*DLL  T 

X\P1=XNT(  IT-3)^DlLT83*(XNPM(3)+XN'PM(1  )  ) -D  EL  T  4  7,  *  X  VP  M  <  ?  ) 
Y\Pl-YNT(lT-3)4.r)rLT83*(  YfjPM(  3>♦Y^PM(1))-JLLT43*Y^PK(?) 
CALL  KG<Tfl«X\Pl,YNPl,XPRM*YPRM,IOUT) 

I P ( I  OUT  .ML  .  C  )GO  TO  1 

XM1  =  X.\T{IT-1)  ♦DELT  3*  (  XPRM^XNPM  (2)  >  ♦DE  L  T  4  3  «  X  MPM  (3  ) 
YMlrYMTdT-D+DELTS^HYPKP-t-YNP*-.  (2>)^DELT43*YMPM(3) 

XMT(IT*i)=XNl 
YMT  <  I  T*  1  )  =YM 
MPOIMTrIT*! 

X  MPM ( 1 ) zxMPM ( r ) 

Xi.'PM  (  2)  =  XMPM  (  7  > 

Y\PM ( 1 ) rYMPM  <  r  ) 

Y;\iPM  (  2)  :YMPM<  3) 

CALL  FG(TM»XM,YM»XMPM(7>,YMPMI3),IDUT) 

I F ( I  OUT .ME . G  )  GO  TO  1 

3  CCNTIMUE 

C  THERE  ./ILL  PL  NPCIMTS  OK  THE  CHARACTERISTIC  LINE 

1  WRITE  (OUTPUTtlOCCOMPOINT 

lOGC  FORMAT ( lx 1 1 ?  ) 

WRlTE<OUTPUT|lCCn((XNT(I)tYNTIl)),I=l*NP;)INT> 

K  J  ]  F  OR  M  A  T  <  b  <  E  X  ,  2  F  1  :  ,  2  )  ) 

RETURN 


METHOD  2  ALGORITHM 
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*****  METHOD  2  ALGORITH!’  ***** 

SUGROUTIM  STREAM  (OUTPUT. X*Y,^ELT) 

SOLVES  F(X»Y.T)  DC/OX  ♦  G(X.Y.T)  DC/OY 
BY  THE  METHOD  OF  CHARACTERISTICS  USING  AN 
RUNGE-KUTTA  FORM  AND  ACCURACY  CHECK. 
INVOLVES  VARIABLE  TIME  STEP  FEATURE. 


DC/DT  =  e 
EMBEDDED  (R.6) 


INTEGER  OUTPUT 

DIMENSION  XNTdDDD.YNT  (ICni  ) 
PARAMETER  TCLA=:‘.!i01 
parameter  TOLBrO.OOOl 
PARAMETER  SMALL  =  1  .E-1  :• 


P ARAMETER 
P ARAMET  ER 
PARAMETER 
I  OUTrC 
WHEN  lOUT 
A  BOUNDARY 
N  P  0  I  N  T  =  1 
INITIALIZE 
X  N  T  (  1  )  =  X 
YNT< 1 )=Y 
XNP  =  X 
YNPrY 
TNP=:  .0 


RETURNED 


FLOW  R'‘GIO!i 


C'^ARACTEP  1ST  If 


CROSSING 


POINT 


ORIGIN 


CHARACTER ISTIC 


I"  ---I 


r.  ■>'.  ■  .  •  . 

I".  "'-I-'-'.-'-l 


I'i! 


PERFORM  RLNGE 

ISTEf  =1 

no  2  ! T  = 1 ,  I  0  C 

r  T  =  DELT/lSTrP 

DT2=LT/2. 

DT4=DT/R. 

DTBrPT/S. 

T  6  =  n  T  /  6  . 

0T27  =  DT  (11, 
DTi3G=DT/3ifc. 
D"’G25=DT/fc2':.. 
DT?3=2.*rT/3. 


RLNGE  -KUTTA 


X  N  =  X  N  P 
YN=YNP 
T  N  =  T  \  P 
KOUNT  =  v 

hithin  each  ^ir, Jor  tine  stea  r)f>lta  T«  Runne-Kutta  is 
used  at  time  steps  delta  t  =  eclta  T/ISTLP. 

DO  3  UT  =  l.IST''p 

CALL  FG(Ti,XN.YN,Fl,Gl.  lOUT) 

I  F  (  I  OUT  .UE  .  0 ) GO  TO  1 

AVAL-.'<N*!)T2*F  1 
YVAL=YN*[ T2*G1 
TyAL=TN*JT2 

CALL  FG(TVALiXVALtYVAL.F2.G2.IOUT> 

I E ( I  OUT .or . n ) GO  to  i 

XV AL  =  XN*r  T  R  ♦  ( F 1 ♦F  2 ) 

YVAl.  =  YN*-T4*<G1^G2) 
tval-tn*dt? 

CALL  FG(TVAL.hVAL»YVAL»F3.G3.10UT) 

I  F  (  TOUT .NE . 0 >00  TO  1 

XVAL=XN4LT*(?.*F3-F2) 

YVAL=YN>DT*(2.*63-62) 


*.  V 

•  li  "  *  "  H 

*>  A  ■  • 


'  s'  s' 

'  s'  s' 

•  .*  s' 


OOO  f-«  IT)  lO  OO  fVi 
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C 
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TVAL=T\>DT 

CALL  FG(TVALtXVAL»YVAL»Fa,GA,ICUT) 

IF (lOUT.NL.C)GC  TO  1 

XVAL  =  XNfDT2  7*(7.*Fl-»10.*P2-*^FA) 

YVAL  =  YN-»DT27*  (7.  ♦Gl+n.* G2  +  G^) 

T VAL=TN*DT23 

CALL  FG(TVALfyVAL»YVALtFG«35*i:UT) 

IFdOUT  .HE.OICO  TO  1 

XVAL  =  XN-*-DT62b*(2«-.*Fl-12‘^.*F2-»5AG.*F3*5A.*F4-3  7R,*F5) 

YVAL  =  YN+PT6  25*(2b.*Gl-12‘:.»G2-»5A6.*G3  +  6A.*G4-3  76  .*Gb) 

TV ALrTN+DT5 

call  FG<TVAL«XVAL»YVAL»Fo*G6,I0UT) 

IF(10UT.NT.3)G0  tq  1 

XN14  =  XN-»CT6*(F1+4.*F3>F<») 

Yrn4  =  YN>DT6*(ri  +  4.*G3+G<») 

XMf,=  XN4-DT3  3fc*(14,*Fl-*.35.*F4-*-lfc2.*F?+125.*F6) 

YM6  =  YM*DT33f.  *<14,*G1*35.*G4  +  162.*G5*125.*6  6> 

TMl=TN*nT 

Peg1n  calculation  for  accuracy  check. 

X  A6=  AhS  <X('J1G) 

Y  AC  =  ABS  <  YN’lfc) 

IF  (XAb.LT.SMALDXAferSMALL 
IF (YAe.LT.SMALL) YA6=CMALL 
CHECKX=A0S(XN16-XN14)/XAF 
CHECK  Y=AHS(YN1G>YM4)/Y  AG 

IF(CH£CKX,LT. TOLA. AND. CHECKY.LT.TOLATGO  TO  5 

Here  if  accuracy  is  insufficient.  Go  back  and  start 
from  the  beginning  of  the  major  time  step  with  a  new 
nelta  t  eoual  to  i/2  the  current  relta  t. 

ISTEP=ISTEP*2 
IF(UTEP.LT.32)GO  TO  4 
WRITE(1,10C2) 

F0RMAT(/lX»3aHMAX  TINE  STEP  TOO  LARGE*  ABORT) 

STOP 

I  F  (CHECK  X  .LT,  TOLL  .AND.CHFCKY  .L  T  .  T  OL  E  )  K  C  UN  T  r  K  G  U  N  T  1 

X\-XM6 

YN=YN16 

TN=TN1 

X:<IT(IT-»1)=XN1G 

YNT(IT-fl)=YNlG 

NP0INT=IT+1 

If  accuracy  criteria  was  consistently  surpassed  throuohout 

the  major  time  step*  double  delta  t. 

lF(K0UNT.F(j.lCTtP)ISTEP-lSTEP/2 

IF( ISTEP.LT.l) ISTEP=1 

XNPrXNl 6 

YNPrYNlG 

TNP  =  TM 

CONT I NUE 

THERE  WILL  BE  NPOINTS  ON  THE  CHARACTERISTIC  LINE 

WRITE(OUTPUTflOCC)NPOINT 

FORMAT (  IX, I  5) 

WRITF(OUTPUT,1001)((XNT(I>,YNT(I)),I=l»NPOlNT) 

FORMAT  (  5<  bX  ,2F1  Cf  .2  )  ) 


•-•ronoo  ooo  o  o  oon*-*  o  n  oono 


FLO*.  LINE  PLOTTFFv 


Paor 


1 


C  »**..  SFLOT  PLOTTING  ROUTINE  •**** 

ROUTINE  FOR  PLOTTING  F  T  R  r M  L  I  N  F  S  .  CREATES  OVERLAYS  WHEN 
UFtO  WITH  APPROPRIATE  SCALE  FACTORS. 

^■NEMONIC  COCES  E  OR  CRREL  FILE  CALLS- 

INTEGER*^  JSHr:AO,J$UPIT»J$RDUR*J$CLOS.J$FRUD»J$RWND 
F’  iRAF  ETER 

1  JJREAD  =  :l,  /*  READ 

1  JSWRIT  =  ;?♦  /♦  WRITE 

1  JSRD.R  =  :3»  /*  READ  AND  WRITE 

1  JJCLOS  =  :R.  /•  CLOSE  file 

1  JSDELF  =  ;5,  /*  DELETE  FILE 

1  JIFRWU  =  :6»  /*  MOVE  FORWARD  GNE  FILE  (MAG  TAPE) 

1  JiRWN'L.  =  :7  /♦  REWIND 

J  I  M  E  N  S I  C  V  X  ( 1  .  0  3  )  ♦  Y  <  1  C  0  ’,  )  ♦  L  A  fc  E  L  C  1  :  ) 

INFILE=5 

CALL  GFILL$(INFILE»JSREAD»»G1  VE  STREAMLINE  INPUT  FILE»,26) 
REAu(INriLE*i:C3)0ELT 
READ  (  INFILL  ;C0>  XI  j  .YM 

000  F0RMAT(lXt4FlD.2) 

INPUT  SCALE.  1  INCH  OF  PAPER  EOUALS  (SCALE)  F E JT , H E T E R S 1 LES t 
PRINT  lOSl 

j;i  FGRMAT{/1X,11PINPUT  SCALD 
KrAD( 1«*)SCALE 

PRINT  1  :  u 

OCA  FORMAT( /lX.31iSYMPOL  EVERY  N  TH  POINT.  GIVE  N) 

PEAU(:«*)INCR 
T INCRxDEL  T* INCR 

INITIATE  PLOT  MODE. 

CALL  PLOTSIv.) 

DETERMINL  FLOW  REGION  DIMENSIONS  IN  INCHES. 

XLONG=( XN-Xl ) /scale 
YLONG=( YM-YI ) /SCALE 

CRAW  OUTLINE  OF  FLOW  REGION. 

CALL  PLOT  (D  -3  ) 

CALL  PLOT (XLONG.C. .2) 
call  PLOT (XL0NG.YL0NG,2  ) 

CALL  PLOT (C.» YLON&.2) 
call  PL0T(C..  :..2) 

READ  IN.  AND  PLOT  THE  STREAMLINES  ONE  AT  A  TIME. 

PEAD(INFILE.i:02.END=l)NPCINT 
CL2  FORMAT ( IX . I F) 

X (NFO INT  +  1 ) =X  1 
X  (NPOINT^D  =SCALF 
Y (NPCINT+1  )=Y1 
Y (NPO  TNT  +  2 ) =SCALE 

?EAD(INFILE,i:C3)((X(I).Y(l)).I=l .NPOTNT) 

PGRMAT('^(5X.2F10.2)) 


1  C  A  3 
C 


i  EKFOFM  PLOTTING  OF  THE  CURRENT  ^LOW  LINE 
CALL  LI Nl (X.Y.NPOlNT.l.  INCR.F) 


>-•000 


FLOW  LINE  PLOTTER 
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<1 


GO  TO  2 


'^UT  REFERENCE  MARKS  AND  ANNOTATION  ON  THE  PLOT. 


E\CODE(18,200C*LAPEL)X1,Y1 
200G  F0RMAT(1H(,F7.1»1H**F7.I*2H)  ) 

SIZE=.16 

yL0C=-.5 

YL0C=-.5 

CALL  SYMBOL (XL0C»YL0Cf SIZE,LAPEL.0.*1P  > 
tNCODE(2C,2C01,LAfiEL)TINCR 
2:01  FDRMATd^HTIMF  INTERVAL  =  *Ff..l) 
XL0C=XL0NG/2. 


CALL  SYMBOL  <XL0C«YL OC *51 2EtLAbELtO.,2r) 
ENCO DEI  1ft *200  0, LABEL) XN,YM 
XL0C=XL0NG-1.75 
YLGC=YL0NG+.5 

CALL  SYMBOL ( XLO C » YLOC ♦$ I 2E *L ADEL ♦ C . ♦ 1  ft  ) 
CALL  FILESi(INFILr,JSCLOS) 

STOP 

END 


